********************************************************************************
******* Miscellaneous results form CBA


************************************************	
clear all
do "$maindir/Code/bootstrapres.do"
use "$maindir/Results/Model_Outputs/cba_data.dta"

*********** Classification analysis - comparing least and most cost-effective homes
sort Household meterend

** usage by month of year
bysort Household end_month : egen tempuse_pre = mean(total_mmbtu) if treated==0 & twoyears_prepost==1
bysort Household end_month : egen monthuse_pre = max(tempuse_pre)
bysort Household end_month : egen tempuse_post = mean(total_mmbtu) if treated==1 & twoyears_prepost==1
bysort Household end_month : egen monthuse_post = max(tempuse_post)

keep Household end_month monthuse_pre monthuse_post ExistingConsumption ProjectedConsumption
duplicates drop
drop if monthuse_pre==. & monthuse_post==.

gen aux = 1
bysort Household : egen tempre = total(aux) if monthuse_pre!=.
bysort Household : egen nmonths_pre = max(tempre)
bysort Household : egen tempost = total(aux) if monthuse_post!=.
bysort Household : egen nmonths_post = max(tempost)

** usage by year
bysort Household : egen mmbtu_pre = total(monthuse_pre)
replace mmbtu_pre = . if nmonths_pre==.
bysort Household : egen mmbtu_post = total(monthuse_post)
replace mmbtu_post = . if nmonths_post==.

keep Household mmbtu_pre mmbtu_post nmonths_pre nmonths_post ExistingConsumption ProjectedConsumption
duplicates drop

gen ratio_pre = (ExistingConsumption/1000000)/(mmbtu_pre)
replace ratio_pre = . if nmonths_pre!=12
gen ratio_post = (ProjectedConsumption/1000000)/(mmbtu_post)
replace ratio_post = . if nmonths_post!=12

foreach x in ratio_pre ratio_post {
	winsor2 `x', trim cuts(0.5 99.5) replace
}

keep Household mmbtu_pre mmbtu_post nmonths_pre nmonths_post ratio_pre ratio_post
duplicates drop

* merging with CBA information based on marginal social costs
merge 1:1 Household using "$maindir/Results/Model_Outputs/cba_scc.dta", keep(3) nogen

* merging state date
merge 1:1 Household using "$maindir/Data/ihwap_state.dta", keep(3) nogen


* merging with information about PRISM engineering model
tempfile ihwapgap
save `ihwapgap'

clear
use "$maindir/Results/Model_Outputs/prism_gas.dta"

keep Household max HDDbest treated
duplicates drop


*** prism-related variables
gen rsq_pre = max if treated==0
bysort Household : egen rsquare_pre = max(rsq_pre)
gen rsq_post = max if treated==1
bysort Household : egen rsquare_post = max(rsq_post)

gen hd_pre = HDDbest if treated==0
bysort Household : egen hdd_pre = max(hd_pre)
gen hd_post = HDDbest if treated==1
bysort Household : egen hdd_post = max(hd_post)

keep Household rsquare_pre rsquare_post hdd_pre hdd_post
duplicates drop

merge 1:1 Household using `ihwapgap', keep(2 3) nogen


********** variables that define low versus high performers
keep if nmonths==12
sort year30_npv
cumul year30_npv, gen(pctile_year30)
replace pctile_year30 = round(pctile_year30, 0.01)
tostring pctile_year30, replace force format(%9.2f)

gen lowperform = 0 if pctile_year30>"0.25" & pctile_year30<="0.75" & pctile_year30!="."
replace lowperform = 1 if pctile_year30<="0.25" & pctile_year30!="."
gen highperform = 0 if pctile_year30>"0.25" & pctile_year30<="0.75" & pctile_year30!="."
replace highperform = 1 if pctile_year30>"0.75" & pctile_year30!="."


* rescaling the NPV variable
gen npv1000 = year30_npv/1000

* indicators for low PRISM r-squared
gen lowrsq_pre = rsquare_pre<0.95 if rsquare_pre!=.
gen lowrsq_post = rsquare_post<0.95 if rsquare_post!=.


***** share of nonzero spending
foreach x in AirCon	AirSeal Attic Baseload ///
			Door Foundation Furnace ///
			General HealSfty WallIns ///
			Window WtHtr {
gen spent`x' = 1 if Real_tot_act`x'>0 & Real_tot_act`x'!=.
replace spent`x' = 0 if Real_tot_act`x'==0
}

gen tab_builddateX = 0 if builddate!=.
replace tab_builddateX = 1 if tab_builddate11==1 | tab_builddate12==1 | tab_builddate13==1

gen liheap = 0 if liheapemergency=="False"
replace liheap = 1 if liheapemergency=="True"

gen pct_BlowerReduc = BlowerReduc/Blower_Pre

gen atticR = Attic_RValue
replace atticR = . if atticR==0

gen hasattic = 0 if Attic_RValue==0
replace hasattic = 1 if Attic_RValue>0 & Attic_RValue!=.

winsor2 sqfeet, replace trim cuts(1 99)


foreach x in female renter white black hispanic nativeamerican otherrace haselderly ///
		hasminor liheap pct_BlowerReduc hasattic nstories tab_builddate1 tab_builddate2 ///
		tab_builddate3 tab_builddate4 tab_builddate5 tab_builddate6 tab_builddate7 tab_builddate8 ///
		tab_builddate9 tab_builddate10 tab_builddateX py09 py10 py11 py12 py13 py14 py15 py16 {
replace `x' = 100*`x'
}

gen projected_savings = 100*engSave
gen realized_savings = 100*realized_savings_mmbtu/counterfactual_mmbtu
gen savings_gap = realized_savings - projected_savings

forval i = 1/200 {
gen realized_savings`i' = 100*realized_savings_mmbtu`i'/counterfactual_mmbtu`i'
gen projected_savings`i' = projected_savings
gen savings_gap`i' = realized_savings`i' - projected_savings`i' 
}

gen ExistingConsumption_mmbtu = ExistingConsumption/1000000
gen ProjectedConsumption_mmbtu = ProjectedConsumption/1000000
 

**** attic spending below and above 1500
gen highattic = 0 if Real_tot_actAttic<1500
replace highattic = 1 if Real_tot_actAttic>=1500 & Real_tot_actAttic!=.

**** furnace replacements (spending above 1800)
gen furnace_replace = 0 if Real_tot_actFurnace<1800
replace furnace_replace = 1 if Real_tot_actFurnace>=1800 & Real_tot_actFurnace!=.


********************************************************************************
********** loop for non-zero spending

*** spent on furnace replacements versus furnace repairs
gen furnspent_replace = Real_tot_actFurnace if furnace_replace==1 & Real_tot_actFurnace>0 & Real_tot_actFurnace!=.
gen furnspent_repair = Real_tot_actFurnace if furnace_replace==0 & Real_tot_actFurnace>0 & Real_tot_actFurnace!=.
label var furnspent_replace "Furnace Replacement"
label var furnspent_repair "Furnace Repair"


label var Real_nonHS_cost "Total, Excluding H&S"
label var Real_tot_actAirCon "Air Conditioning"
label var Real_tot_actAirSeal "Air Sealing"
label var Real_tot_actAttic "Attic Insulation"
label var Real_tot_actBaseload "Baseload"
label var Real_tot_actDoor "Door"
label var Real_tot_actFoundation "Foundation"
label var Real_tot_actFurnace "Furnace"
label var Real_tot_actGeneral "General"
label var Real_tot_actHealSfty "Health and Safety"
label var Real_tot_actWallIns "Wall Insulation"
label var Real_tot_actWindow "Window"
label var Real_tot_actWtHtr "Water Heater"

foreach x in Real_nonHS_cost Real_tot_actAirCon	Real_tot_actAirSeal Real_tot_actBaseload ///
			Real_tot_actDoor Real_tot_actFoundation furnspent_repair furnspent_replace ///
			Real_tot_actGeneral Real_tot_actHealSfty ///
			Real_tot_actWindow Real_tot_actWtHtr ///
			Real_tot_actAttic Real_tot_actWallIns {
	sum `x' if highperform==0 & `x'!=0
	sca coef1 = `r(mean)'
	sum `x' if lowperform==1 & `x'!=0
	sca coef2 = `r(mean)'
	sum `x' if highperform==1 & `x'!=0
	sca coef3 = `r(mean)'
	reg `x' i.lowperform if `x'!=0
	matrix mat1 = e(b)
	sca coef4 = mat1[1,2]
	reg `x' i.highperform if `x'!=0
	sca nobs = e(N)
	matrix mat2 = e(b)
	sca coef5 = mat2[1,2]
	matrix ttestcoefs = (coef1, coef2, coef3, coef4, coef5)
	mat colnames ttestcoefs="A" "B" "C" "D" "E"

	matrix input btstcoefs = (., ., ., ., .)
	forval i = 1/200 {
	qui : sum `x' [fweight=weight_boots`i'] if highperform==0 & `x'!=0
	sca coef1 = `r(mean)'
	qui : sum `x' [fweight=weight_boots`i'] if lowperform==1 & `x'!=0
	sca coef2 = `r(mean)'
	qui : sum `x' [fweight=weight_boots`i'] if highperform==1 & `x'!=0
	sca coef3 = `r(mean)'
	qui : reg `x' i.lowperform [fweight=weight_boots`i'] if `x'!=0
	matrix mat1 = e(b)
	sca coef4 = mat1[1,2]
	qui : reg `x' i.highperform [fweight=weight_boots`i'] if `x'!=0
	matrix mat2 = e(b)
	sca coef5 = mat2[1,2]
	matrix tempmat = (0, 0, 0, coef4, coef5)
	matrix btstcoefs = (btstcoefs \ tempmat)
	}
	mata : st_matrix("ttestsvar", mm_colvar(st_matrix("btstcoefs")))
	matrix ttestsvar = diag(ttestsvar)
	mat colnames ttestsvar="A" "B" "C" "D" "E"
	mat rownames ttestsvar="A" "B" "C" "D" "E"

	bootstrapres ttestcoefs ttestsvar nobs
	est store `x'
}

esttab Real_nonHS_cost Real_tot_actAirCon	Real_tot_actAirSeal Real_tot_actBaseload ///
			Real_tot_actDoor Real_tot_actFoundation furnspent_repair furnspent_replace  ///
			Real_tot_actGeneral Real_tot_actHealSfty ///
			Real_tot_actWindow Real_tot_actWtHtr ///
			Real_tot_actAttic Real_tot_actWallIns , se nostar
matrix C = r(coefs)
eststo clear
local rnames : rownames C
local models : coleq C
local models : list uniq models
local i 0
foreach name of local rnames {
    local ++i
    local j 0
    capture matrix drop b
    capture matrix drop se
    foreach model of local models {
        local ++j
        matrix tmp = C[`i', 2*`j'-1]
        if tmp[1,1]<. {
			local varlab : variable label `model'
            matrix colnames tmp = "`varlab'"
            matrix b = nullmat(b), tmp
            matrix tmp[1,1] = C[`i', 2*`j']
            matrix se = nullmat(se), tmp
        }
    }
    ereturn post b
    quietly estadd matrix se
    eststo `name'
}
esttab B A C D E using "$maindir/Results/Tables/ttests_nonzerospending_year30.tex", ///
	replace se b(3) se(3) noobs nonotes ///
	mlabels("Bottom 25% Homes" "Interquartile Homes" "Top 25% Homes" "Diff. (1)-(2)" "Diff. (3)-(2)")
	
	
********************************************************************************
**** loop over covariates related to energy savings 
label var ratio_pre "Ratio Modeled/Actual Usage Pre"
label var ratio_post "Ratio Modeled/Actual Usage Post"
label var realized_savings "Realized Savings (%)"
label var projected_savings "Projected Savings (%)"

forval i = 1/200 {
	gen ratio_pre`i' = ratio_pre
	gen ratio_post`i' = ratio_post
}

foreach x in ratio_pre ratio_post realized_savings projected_savings {
	sum `x' if highperform==0
	sca coef1 = `r(mean)'
	sum `x' if lowperform==1
	sca coef2 = `r(mean)'
	sum `x' if highperform==1
	sca coef3 = `r(mean)'
	reg `x' i.lowperform
	matrix mat1 = e(b)
	sca coef4 = mat1[1,2]
	reg `x' i.highperform
	sca nobs = e(N)
	matrix mat2 = e(b)
	sca coef5 = mat2[1,2]
	matrix ttestcoefs = (coef1, coef2, coef3, coef4, coef5)
	mat colnames ttestcoefs="A" "B" "C" "D" "E"

	matrix input btstcoefs = (., ., ., ., .)
	forval i = 1/200 {
	qui : sum `x'`i' [fweight=weight_boots`i'] if highperform==0
	sca coef1 = `r(mean)'
	qui : sum `x'`i' [fweight=weight_boots`i'] if lowperform==1
	sca coef2 = `r(mean)'
	qui : sum `x'`i' [fweight=weight_boots`i'] if highperform==1
	sca coef3 = `r(mean)'
	qui : reg `x'`i' i.lowperform [fweight=weight_boots`i']
	matrix mat1 = e(b)
	sca coef4 = mat1[1,2]
	qui : reg `x'`i' i.highperform [fweight=weight_boots`i']
	matrix mat2 = e(b)
	sca coef5 = mat2[1,2]
	matrix tempmat = (0, 0, 0, coef4, coef5)
	matrix btstcoefs = (btstcoefs \ tempmat)
	}
	mata : st_matrix("ttestsvar", mm_colvar(st_matrix("btstcoefs")))
	matrix ttestsvar = diag(ttestsvar)
	mat colnames ttestsvar="A" "B" "C" "D" "E"
	mat rownames ttestsvar="A" "B" "C" "D" "E"

	bootstrapres ttestcoefs ttestsvar nobs
	est store `x'
}

esttab ratio_pre ratio_post realized_savings projected_savings, se nostar
matrix C = r(coefs)
eststo clear
local rnames : rownames C
local models : coleq C
local models : list uniq models
local i 0
foreach name of local rnames {
    local ++i
    local j 0
    capture matrix drop b
    capture matrix drop se
    foreach model of local models {
        local ++j
        matrix tmp = C[`i', 2*`j'-1]
        if tmp[1,1]<. {
			local varlab : variable label `model'
            matrix colnames tmp = "`varlab'"
            matrix b = nullmat(b), tmp
            matrix tmp[1,1] = C[`i', 2*`j']
            matrix se = nullmat(se), tmp
        }
    }
    ereturn post b
    quietly estadd matrix se
    eststo `name'
}
esttab B A C D E using "$maindir/Results/Tables/ttests_save_year30.tex", ///
	replace se b(3) se(3) noobs nonotes ///
	mlabels("Bottom 25% Homes" "Interquartile Homes" "Top 25% Homes" "Diff. (1)-(2)" "Diff. (3)-(2)")	


	
***************************************************************
********** histogram of gap by benefits
	
gen pct_gap25 = savings_gap if lowperform==1
gen pct_gap75 = savings_gap if highperform==1
gen pct_gap50 = savings_gap if highperform==0

gen gap_bins = .
forval i = -10(10)60 {
	replace gap_bins = `i' if savings_gap>`i'-10 & savings_gap<=`i'
}
replace gap_bins = -20 if savings_gap<-20
replace gap_bins = 70 if savings_gap>=60 & savings_gap!=.

gen gapbinlab = "<-20%" if gap_bins==-20
replace gapbinlab = "[-20% -10%)" if gap_bins==-10
replace gapbinlab = "[-10% 0%)" if gap_bins==0
replace gapbinlab = "[0% 10%)" if gap_bins==10
replace gapbinlab = "[10% 20%)" if gap_bins==20
replace gapbinlab = "[20% 30%)" if gap_bins==30
replace gapbinlab = "[30% 40%)" if gap_bins==40
replace gapbinlab = "[40% 50%)" if gap_bins==50
replace gapbinlab = "[50% 60%)" if gap_bins==60
replace gapbinlab = ">=60%" if gap_bins==70

labmask gap_bins, values(gapbinlab)

sum savings_gap

graph bar (count) pct_gap25 pct_gap50 pct_gap75 ///
		, over(gap_bins, label(angle(45) labsize(small))) stack graphregion(color(white)) bgcolor(white) ///
		legend(style(column) title("Home Cost-Effectiveness:", size(small)) ring(0) position(2) bmargin(medium) ///
		order(1 "Bottom Quartile" 2 "Interquartile" 3 "Top Quartile") ///
		symysize(5)  symxsize(5)) bar(1, color(maroon)) bar(2, color(navy)) bar(3, color(forest_green)) ///
		ytitle("Number of Homes") b1title("Performance Wedge (%)" "(realized - projected savings)")
graph export "$maindir/Results/Figures/gap_histogram_year30.png", replace width(2500)

